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A FINITE ELEMENT METHOD FOR HIGH-CONTRAST INTERFACE 
PROBLEMS WITH ERROR ESTIMATES INDEPENDENT OF CONTRAST 

JOHNNY GUZMANS MANUEL A. SANCHEZS AND MARCUS SARKIS^ 


Abstract. We define a new finite element method for a steady state elliptic problem with 
discontinuous diffusion coefficients where the meshes are not aligned with the interface. We 
prove optimal error estimates in the norm and weighted semi-norm independent of the 
contrast between the coefficients. Numerical experiments validating our theoretical findings are 
provided. 


1. Introduction 

In this article we develop a finite element method for a steady state interface problem. We 
pay particular attention to high-contrast problems, proving optimal error estimates independent 
of the contrast of the discontinuous constant coefficients for the numerical method. 

Let If C be a polygonal domain with an immersed interface F such that = 17 U 17 
with 17“ nl7+ = 0, and F = 17 n 17"*". We assume that F does not intersect 917, enclosing either 
17“ or 17"*“. Our numerical method will approximate a solution of the problem below. 

(Fla) —in 17^, 

(1.1b) u = 0 on 917, 

(1.1c) [u] =0 on F, 

(l.ld) [pDnu] =0 on F. 

The jumps across the interface F are defined as 

[pDnu] = p~D^-u~ + p~^D^+u~^ = p~Vu~ ■ n~ + p'^Vu'^ ■ n^, [u] = u~^ — u~, 

where = n|Q± and is the unit outward normal to 17^. We furthermore assume that 

p'*' > p“ > 0 are constants and that the interface F is a closed, simple and regular curve with 
an arc-length parameterization X. 

There has been a recent surge in the development of finite element methods for interface prob¬ 
lems. See for instance [2111210 [3 0 [la 0 iig [la iig [Ill |22i 13 ia E [g to name a few. Among 
the articles where the discretization is based on meshes not aligned with the interface, most of 
the methods focus on low contrast problems and only a few address the high contrast problems 
(p'*“/p“ ;g> 1). For example, Burman et al. [7] introduced an unfitted Nitsche’s method with 
averages and stabilization techniques for arbitrarily high-contrast problems, presenting bounds 
for the condition number of the stiffness matrix, although a rigorous error analysis was not 
given in that paper. Another example, is given by Chu et al. [8] that uses multiscale techniques 
to build basis functions, an approach that seems well suited for high curvature problems (e.g. 
inclusions completely contained in a triangle). However, in regions where the curvature of the 
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interface is small it appears that their main a priori estimate degenerates (see Theorem 3.9 in 
m), forcing them to refine the mesh on those regions in order to be aligned with the interface. 
In our approach we do not need mesh refinements to make the triangulations aligned with the 
interface, however, we do not address high curvature problems. 

In order to put our contribution in context, let us explain two popular finite element ap¬ 
proaches for problem (1.1). The first approach is to double the degrees of freedom on triangles 
that intersect the interface and then add penalty terms to weakly enforce the continuity across 
the interface, see for example [7]. Burman et al. [7] demonstrated that in addition to penalizing 
the jumps across T it is necessary to add a flux stabilization term. This method is the so-called 
stabilized unfitted Nitsche’s method. The stabilization term penalizes the jumps of the gradient 
on edges that belong to triangles that intersect the interface. As Burman et al. [7] showed, 
in order to obtain a method that is robust with respect to diffusion contrast and robust with 
respect to the way T cuts triangles, this type of penalization is necessary. The second common 
approach, and the one we focus in this paper, is to define local piecewise polynomial finite el¬ 
ement spaces on triangles that intersect the interface T (see for instance m na m El Eg IT]). 
The basis functions are constructed by having them satisfy the continuity of the solution and 
the continuity of the flux strongly across T. Unlike the unfitted Nitsche’s method that weakly 
imposes the interface conditions, in this approach the flux conservation and the continuity of the 
solution are enforced strongly, for example at certain points on T, without requiring stabilization 
terms on T. This is an important and distinguishable feature of the Immersed Interface Method 
using a Finite Element formulation (Immersed Finite Element Methods, see [20) 1. These basis 
functions are defined locally on each triangle, and therefore they are naturally discontinuous 
across edges of the triangulation. Namely, Adjerid et al. in [1] proposed to penalize jumps of 
the trial functions across the edges. Similarly, Lin et al. in [2T] added similar penalty terms and 
proved optimal error estimates. However, in their analysis they do not consider high-contrast 
problems. 

In this paper we follow this approach, defining local basis functions that are piecewise poly¬ 
nomials on each side of triangles that are cut by T. However, an additional stabilization term 
is added as compared to the methods of Adjerid et al. [I] and Lin et al. [2T], allowing us to 
prove error estimates that are independent of the contrast jp~. By stabilization term we 
refer to a penalization of the jumps of the normal derivatives of the approximation across the 
edges that belong to triangles intersecting T. This idea was used before by Burman et al. [7], 
however here we use different stabilization parameters (and of course different basis functions). 
Roughly speaking, the reason this flux stabilization is important for high contrast problems, is 
that one does not want to move estimates from n"*" to 0“ because p'^ could be much larger than 
p~. However, triangles that are cut by T might have a thin part in n"*" and therefore inverse 
estimates might be affected. By adding the jumps of derivatives we can transfer the estimates 
to a neighboring triangle that will have a larger portion in 

In order to prove error estimates independent of the contrast p'^ jp~ we assume tP' regularity 
of u on both and 0“. To be more precise, the error estimate presented in this paper for the 
energy norm (see (4.1)) is 


( 1 . 2 ) 


\u-Uh\\v {\\Du\\l2^q-) + \\D‘^u\\l2^q-)) 

+\/^(||T*u||i2(f^+) -I- ||L>^tt||L2(Q+))| 
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Consequently, assuming that 0, is convex and using regularity estimates (see 0) we will have 
the result 


(1.3) 


u — UhWv 


< C 



ll/llL2(n)- 


In addition, using a duality argument we prove the estimate \\u — Uh\\i 2 (^Q^ < C{h‘^/p~)\\f\\]^ 2 (^Qy 
The constants in the estimates depend on the geometry, including the curvature of T. 

The outline of the paper is as follows. We formulate the discrete problem as finding Uh G 14 
such that ah{uh,Vh) = {f,Vh), Vu/i G I 4 . The discrete space I 4 is introduced in Section 2.1 
and the bilinear form a/j(-, •) in Section 2.2 In Section fundamental resnlts on element¬ 
wise weighted and norm approximation for the space I 4 are established. Coercivity 
and continuity of the bilinear form ah are studied in Section 4.1 and Section 4.2, respectively. 
In particular for the continuity, we note that the use of an augmented norm is necessary for 
the analysis due to the presence of the penalty terms involving flux jumps. In Section the 
bound (1.2) is established by estimating the approximation error and the consistency error 
across T and across elements near T. The error estimate in the norm is also proved in this 
section. Section is devoted to present extensions of the method to three dimensions, discnss 
related methods, and state some concluding remarks. Finally, in Section we provide nnmerical 
experiments corroborating our theoretical findings. An Appendix containing technical proofs 
and a computational consideration is also included. 


2. The finite element method 


2.1. Notation and local finite element space. In this section we present a hnite element 
method for problem (1.1) nsing piecewise linear polynomials. 

We next develop notation. Let Th, 0 < h < 1 be an admissible family of triangulations of 
O (conforming), with fl = DTeTh'^ the elements T are mutually disjoint. Let hx denote 
the diameter of the element T and h = maxj’ Ht- We let be the set of all edges of the 
triangulation. We adopt the convention that edges e, elements T, sub-edges := e H 
sub-elements := T n and sub-regions are open sets, and we nse the over-line symbol 
to refer to their closure. 

Let denotes the set of triangles T ^ Th such that T intersects T. We let be the set of 
all the three edges of triangles in . Since T is we have that ||X"|4oo < cx) and we dehne 
the maximum curvature k := Our analysis, will be valid when h is sufficiently small. 

To make this precise, we will use the concept of a tabular neighborhood whose existence is a 
standard resnlt in differential geometry; see [H] Section 2-7, Proposition 3. 


Lemma 1. (Existence of r-tubular neighborhood) Let T be a regular, simple, curve. For 
every x G T consider the line segment Nx{r) of length 2r centered at x and perpendicular to T 
at X. Define the tubular neighborhood of radius r 0 /T by Tub{r) = UxerLfxir). Then, there 
exists r > 0 such that for any two points x,y £ T, x ^ y, the line segments Nx{r) and Ny{r) are 
disjoint. 


From now on we will work under the following assumptions. 

Assumption. Given Lemma\^we make the following assumptions on the triangulation: 

(1) We assume that the triangulation is shape-regular, see [3]. 

(2) We assume h < r/2 where r is the radius of the tubular neighborhood o/F. 

(3) The interface intersects the boundary of an element at most twice and at different edges. 









4 


JOHNNY GUZMANi, MANUEL A. SANCHEZ^, AND MARCUS SARKIS^ 


It is well known that r < 1 /k, and hence by our assumption h < 1/(2k). The radius r also 
bounds from below how close the curve T comes from self-intersecting (e.g. consider a dumbbell 
with a thin middle section). We make use of these assumptions to prove some of the technical 
lemmas (see Lemmas §[Tg[g. 

In addition, we dehne an element patch ujt of a triangle T, its restriction to and its 
intersection with T by 

ut ■= Int{ 1^ K : K nT ^ 0}, := Wy := wt H T, 

K&Th 

where “Int” denotes the interior of the set. 

The introduction of notation for the patch will be relevant in the proof of the interpolation 
error further forward. We first need to build our local hnite element space (on each element). 
To this end, we let xq be the midpoint with respect to the arc-length on the curve segment 
Tr := T n T. We note that the midpoint choice is a preference of the authors, the proofs below 
hold for any xq G Tp. Let Lt be the line segment inside T which is tangent to T at xq. We 
define by the unit tangent vector to T by a 90° clockwise rotation of n^. Figure [^illustrates 
the definitions and notations introduced above. 


Figure 1. Illustration 
of our notation on 
an element T G . 



In order to define our finite element space, we will need the following lemma. 
Lemma 2. Consider the operator T ; P^(T+) — )• P^(T“) defined by 


( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 


T(x)(xo) 

(^t+^(^))(a^o) 

P~iD^+T{v)){xo) 


x(xo), 

{D^+v){xo), 

P^iDr,+ v){xo), 


where rig = n=*=(xo) and tg = t^(xo). Then, T is well defined. 

Note that the exact solution satisfies the transmission conditions (2.1), (2.2) and (2.3) on 
all points over F. 

Next, given T G and for each v G P^(r'’') we can consider the unique corresponding 
function 


Giv) = 


V, 

T(i;), 


on T+, 
on T~. 
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Let span {vi,V 2 -,v^] be a basis for P^(T) restricted to T~^. Then we define the local finite 
element space 


(2.4) 


S'i(r) = [span {G(?;i),G(i;2),G(i; 3)}, if T e , 


¥\T), iiTeThW^. 


An explicit construction of basis functions of the local space S^{T) is given in Appendix]^ The 
global finite element space is defined by 

Vh '■= {v '■ v\t ^ VT G 7h, u is continuous across all edges in . 

2.2. Finite element method. We begin this section by introducing some standard discontin¬ 
uous finite element notation for jumps and averages. 

For a piecewise smooth function v with support on 7/i, we define its average and jump across 
an interior edge e G S^\dQ, shared by elements Ti and T 2 , as 

i^Iti + I^|t 2 


w = 


H = ulrjni -h v\T2n2, 


where rii and n 2 are the outward pointing unit normal vectors to Ti and T 2 , respectively. 
Similarly, if r is a vector-valued function, piecewise smooth on Th, its average and normal jump 
across an interior edge e are defined as 




tIti + t\t2 


M = 't\ti ■ ni -h t\t 2 ■ rL2. 


Next we introduce the finite element approximation to problem (1.1). First, we define the 
space V as the union of the broken Sobolev spaces 

(2.5) V := U where hI{Q^) = {v : v\t± G for all T G %}■ 

We can then define the bilinear form : F x F —?■ M and the linear functional (/,•): F —?• M by 

(2.6) ah{w,v) := j pVh'w-VhV-'^ j ({pV/iu} • [u;]]• [u]]) 

u —p t/ 6 


+ E 
+ E 


7 






M + 


7 


P M 




{f,v) := 

T&Th 


p [V/^u]] [V/iU;]] \e 

rv+ f 


pHVhvj 


T- 


lT+ 




for a penalty parameter 7 > 0. The discrete gradient operator Vh is piecewise defined on 
for an element T G 7/i by 

^ hv\T± = Dv = Vu. 

Finally, the finite element approximation solves: Find Uh € Vh such that 
(2.7) ah{uh,v) = {f,v), for all u G F/j. 


Note that in (2.6) we not only penalize the jumps of the function but also the normal jumps 
of the first derivatives across edges. This will allow us to prove coerciveness and a priori error 
estimates independent of the contrast of the coefficients and also independent of how small T~^ 
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Figure 2. Illustration of basis functions on triangle in Figure with p'^ = 100 
and p~ = 1. 


or T might be. Finally, we like to stress that in (2.6) only the normal derivative jumps are 
penalized, not the tangential derivative jumps. 


3. Local Approximation 


In this section we show that our finite element space has optimal local approximation prop¬ 
erties. Henceforth we will keep track, as much as possible, on how constants depend on the 
maximum curvature k, and in general in the geometry of sub-domains. When we use standard 
results (e.g. regularity and extensions lemmas) we just state that the constants depend on the 
geometry without explicitly saying how they scale with quantities such as curvature and the 
radius of the tubular neighborhood r. In order to define an interpolation operator onto we 
first state an extension result. Consider G (with = 0 on Then, there 

exists a constant C and extensions G H^{Q) with the following properties 

(3.1) Rg = in 

(3-2) 11111,2 (n) + \\D‘^u'^\\L'2(n) < C{\\Du^\\j^2(^q±-^ + \\D‘^u^\\L2(n±))- 

This result follows from Theorem 7.25 in m and applying Poincare’s inequalities. Considering 
that we will only require the extensions to be defined on patches of elements on , in fact, 
we only need the following more local bound which could lead to a better geometric constant: 

(3.3) ||T)ug||E2(f2±|j'r«b(2/i)) + \\D‘^'^%\\L'^{n±\jTub(2h)) < Ce{\\Du^\\l2(q^±'^ -[- \\D‘^u'^\\i2(^^±)). 

For the sake of simplicity, we do not prove here how Ce depends on the geometric constants 
(e.g. K and r). 


Definition 3.1. Let G LI^(H^). For each T G we define Itu G S^{T). In fact, we define 
Itu on all lot 


Itu = 


Il^u, in Lo 
Irp ri. 


m LO 


where are defined satisfying the following conditions 


+ 

T 

r> 


(3.4) 


{i:^u){xo) 

< (D t+I t u)ixo) 

p-{D +1-u){xo) 

\ U 


{Jtu%){xo) 

iIIt+('lTUE)){xo) 

P~iIIn+’lTUE)ixo) 


{i:}:u){xo) 

(D.+/±r)(xo) 

P'^iIIn+lT'^)i^0)- 
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Here Jt is the L?' projection operator onto (note that ojt C Tub{r) since we are assuming 

2h < r). 

If T does not belong to , we consider the sets = U{T £ Th ■ T G and we define 

Itu = IsxU^ 1 for T C where IszU^ is the Scott-Zhang interpolant of defined on 
However, we need to modify the interpolant at every vertex x (if one exists) that is an endpoint 
of an edge e satisfying e C n (e.g. e C T). In this case, we define IszU^{x) to be the 
average of on the edge e. Note that such an x might correspond to two edges. Either edge will 
work for the definition of IszU^{x)- This will give IszW^ = IszU~ for every edge e C H 
since u'^ = u~ on such an edge. 

Consequently, we define the interpolation operator Ih onto the finite element space Vh as the 
restriction of the local interpolation operator It, i.e. 


(3.5) Ihu\T = It[u)\t, for all T E 7)i. 

The local interpolant Itu was constructed so that Lemma holds which in turn is the main 

The interpolator Itu was designed so that if we multi- 
and /?■' 


tool to prove the crucial Lemma [TO 
ply the whole inequalities 

(after using p~ < /)+) will appear 


and (3.7) by p 


respectively, only terms of the form 


Since the proof of Lemma is quite involved, we first prove a local energy stability of the 
interpolant whos proof is much easier and that motivates the definition of Ih- More specifically, 
we prove the following lemma. 


Lemma 3. It holds, 
where C is independent of p^. 


Proof. Clearly we have 

Using that D^+{P^u) and are constant and using the definition of the interpolant 

(3.4) we have 

p+||V(/rn)|li2(^+) < <^(P+Pt+(Jr4)llL2(a;+) +^■|l^r^+('^TI^E)llL2(^+))• 

We trivially have 

p+||V(lTn)||^.(^+) < C(p+||V(Jt4)||l 2(,,^) +p-||V(JTU^)||i2(,,^)). 

Define and note that Jtc^ = c^. We then see 

P+||V(/tu)||^,(^+) < C{p+\\V{JTiu+ - c+))||i2(,,^) + p-\\V{Mu-^ - c-))||i2(,,^)). 

Using an inverse estimate, the stability of the projection Jt, and Poincare’s inequality we 
get 

P+||V(/Tl^)|li2(^+) < C'(/j’^||VM^||i2(^,^) +p“||VUg||i2(^,^)). 

In a similar way we can prove the estimate for ||V(/^tt)||^ 2 (^j-) where we only need to use 

P~ < P^- 

□ 


Lemma 4. Consider T E . Then, for every v E P^(a;T) ihe following hound holds 


/i^||DJu||i2(^^) < C (^hT\v{xo)\ + h‘^\D^+v{xo)\ + h‘^\D^+v 


3 = 0 , 1 . 
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Proof. Using the Cauchy-Schwarz inequality one has 

Using the fact that u is a linear function, the lemma follows. □ 

We next prove a fundamental result of this paper, a local approximation property on the 
space S^{T). 

Lemma 5. Let G satisfying the interface conditions [u] = 0 and [pDnu] = 0. Then, 

for any T G Tj^ and j = 0, 1, the following bounds hold: 

(3.6) 

< C{l + K)hl (II^^^BllL2(a;r) + IIII+ W^'^'^eWl^m) 

and 

(3.7) h^j.\\D^{u'^ — L^u)\\^2(^^+'f ^ C {1 + hi)hrp(^\\Du'^\\];^2(^i^.j.-^ + \\D 

Proof. First note that by adding and subtracting Jtu^ to u^—I^u, using the triangle inequality, 
and by means of the following well known approximation property for (p G H‘^{ujt) 

(3.8) /irll^^(<7>-^T(^)||L 2 (a;j,) < C'^tII^V||l 2 (^,j,) forj = 0, 1, 
the proof of the lemma reduces to estimate 

hlp\\D^w^W where := Jtu^ — I^u, for G P^(a;T)- 

According to the definition of It (3.4), and denoting = n=*=(xo) and tg = t^(xo), we have 

wf{xo) = JTuf.{xo) - Jtu%{xq), 

{D^+wf){xo) = {D^+JTuff){xo)- {D^+Jtu^){xo), 


{D„+w^){xo) = 0 , 


and 


w^{xo) = 0 , 

{D^+wi^){xo) = 0 , 

P'^iDn+'^T)i^o) = P'^{D^+JtU%){xo) - p-{D^+JTUf.){xQ). 

Then, using Lemma and the values derived for we have 

(3.9) ^tII^'^^tIIl2(<.j-) ^ c{hT\wf{xQ)\+h^\D^+wf{xo)\^, 
and 

(3.10) ^tII'^'^^tIIl2(^+) < Chrp\D^+w'^{xQ)\. 

We proceed by bounding the two terms in (3.9) and the term in (|3.10 ), separately. For the first 

- — -- 

Orp, 


term in (3.9), we use that u is continuous on the interface (1.1c), in particular on wS, and we 


apply the triangle inequality 

\wf{xo)\ = \{JTUf; - JTU^)ixo)\ <\\JtUe- 
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A.l 


By means 


Note that cHt < < Chx, where the second ineqnality follows from Lemma 

of a Sobolev inequality in one dimension, we have 

\\JtUe - Ue\\l°°{u>^) < C (^^-j=\\JtUe - «EllL2(^r) + ^At||-D( - '“e)IIl2K)^ , 

consequently, using a trace inequality we obtain 

W-^TUe ~ '^e\\l°°{iJ^) — ^ iP'T ~ W^^’^TUe ~ Ue)\\l‘^(ujj.)) 

Hence, using the approximation property ( |3.8| ), we get 

\\JtUe - Ue\\l^{u,^) < C'/it||-D^«eIIl2(^,j,). 

Analogously, we can show that 

\\JtU% - uy\E^{uj^) ^ ChT\\D‘^U^\\L2(^^^y 


Therefore, we have the bound for the first term in (3.9) 

(3.11) /ir|rcy(xo)| < C/if. (||L»\^||i 2 (^j,) + ||L>\+||i 2 (^,^)) . 

We now turn to the second term in (3.9). We use the fact that D^+w^ is constant on to 
obtain 

l-^t+^r(®o)| = I^^tI < Ch^ ^ W^i+^'^tUe — JTU%)\\E2i^^^Ty 

Using the identities , denoting t+ = t^{x) and n+ = n+(x), for x G wf) 

(3.12) D^+u% = ■t-^)Dt+u% + {t+■n+)D^+u%, 

(3.13) D^+Ue = D^+Ue, 


(3.14) 
we have 


Dr,+ut = ^D„+u 


P 




llL2(^r) = \\D^+{JtUe -Ue + U+- JtU%) + (1 - ^)(*0 • 

< \\D[JtUe - Ue)\\l^(uj^) + \\D{JtU% - U%)\\e2(^V) 

+ (i-^)II(*o+-^+)^«eIIl2 K). 

For the first two terms in the previous bound we use a trace inequality to obtain 

\\D{Jtu% - u|)||i2(^r) < C (^-^\\D{Jtu% - u%)\\l2(^^) + ^/hT\\D‘^{JTU% - w|)||l2 (^^)^ , 

and for the third term we use that \tQ ■ n~\ < ||W"||j;,oo|a;f(| = 0{KhT) on cjf., p~ < , and a 

trace inequality to obtain 


(1 - ^)ll(*o • '<^'^)D'^e\\l^{u2^^) ^ CK,{yh^\\DuE\\L^{i^^) + H^^'^WD'^UeWl^ujt)) ■ 

Hence, applying approximation property ( |3.8[ ) we obtain 

(3.15) /if.|L>^+u;^(xo)| < ^(1 + K)/if. (||L»u“||^2(^,^) + ||i2(^,^)) + ||L»\+||i2(^j,)) . 
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lemma, inequality (3.6). 


If we combine the inequalities (3.11) and (3.15) with (3.9), we arrive at the hrst result of the 


Now we estimate the term in (3.10). We need to bound 

\D^+w+{xo)\ = < Ch~^^‘^\\D^+{JTU+ - ^JTUE)\\L 2 (u^y 


Similarly to the previous bound, we note that 


^n+^E = (^0 • t^)Dt+UE + (rio • n+)D^+u 


E’ 


and then using (3.13) and (3.14) we obtain 


\D 


+ ^TllL2(tjr) — U'^ tt^)) + (1 ){nQ ■ t^)D^+U~^\\^2(^^Ty 


p 


p-r ^ ^ , p-. 

The remaining of the proof is similar as above and we obtain 

^Tl-^n+'“^T(®o)| < C{1 + K)h\' + 1111^,2) + \\D^U%\\^2 


If we combine this inequality with (3.10), we arrive at our second result (3.7) 


□ 


4. COERCIVITY AND CONTINUITY OF BILINEAR FORM 
The aim of this section is to prove coercivity and continuity of the bilinear form ah{-,-) defined 


in (2.6). 


4.1. Coercivity of Bilinear Form. We define the following energy norm || • ||u : V 

(4.1) ii»ii?.. = ll^/pv^^,|lL|„, + ^ ('t((-iihiiL,.-, + 






In order to prove coercivity we will need the following lemma. 

Lemma 6. Let e = Int(9ri n 5 T 2 ) G £\, for Ti,T 2 G T^. Then, there exists a constant 9 > 0 
such that 

< 0max|r,i^|. 
i=l,2 * 


The constant 9 depends on the shape regularity of the triangulation. 
Proof. See Appendix |A.1[ 


□ 


Lemma 7. (Coercivity) If ^ is large enough (depending on shape regularity of triangulation) 
there exists a constant c > 0, independent of h, p~ and p'^, such that 


c||u||y < ah{v,v), for ally £Vh. 


(4.2) 
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Proof. Let v £ Vh, then 

ahiv,v) = |IVpVftu||^ 2 (f^) - 2 ^ [{p'^hvj-M 

J e 




IIHIli 2 (e-) + j^/5’^||HllL2(e+)'j 

^ 

+ E (p“|e-|lllV,nl||i 2 (,-)+p+|e+|||[V,nl||i 2 (,+)). 




To prove the lemma, it is enough to bound the non-symmetric term 

/ {p^hv} • H = / {P^hv} • H + / {P^hv} ■ H- 
Je Je- Je+ 

Let Ti,T 2 G 7ft be such that e = Int(9ri n 8 X 2 ) and set Te = Ti U T 2 . Without loss of 

generality assume |T 2 ~| = max|T~|. Then, according to Lemma 6 

i=l,2 \J 

(4.3) \e-\‘^ <e\Tf-\. 

Then, we see that 

j^JpVv} -Ivj = 

= P~iDn 2 v) {v\t 2 + Dn 2 v){v\T^ -v\t 2 ) 

= [p~{Dn2v) iv\T2 -'f^|Ti) + ^([[Vfti;]])(i;|Ti -v\t2)- 


Using the fact that n is a linear function on T 2 and using (4.3) we have 

1 


(4.4) 


|2 y_p (Dr^^v) {v\t 2 - v\t^)\ < C\\p - Dn ||^2( T -) j ^ rjYy ^|| MllL 2( e -)> 
where C depends on 9. Therefore, we have 

|2 [p~{Dr,^v) iv\T 2 -v\Ti)\ < e/9"||Vfti;||5^2^.^_^ + 


,„2 


e e 


'lllL 2 (e-)’ 


for any e > 0. Furthermore, we have 

I P~{r^hvj) (uIti - w|t 2)| < ep"|e"|||[[Vftn]]||^2(e-) + 
Collecting the last two estimates gives 

|2 [_{pVhv} • HI < ep-{\e-\\\lVhvj\\l2^^.) + ||Vftr;|| 2 ,^^_p + 


e e 


'^llli 2 (e-)- 


+ i)p- 


e e 


^lllL 2 (e-)- 


Likewise, we can bound the integral over e"*" to get a combined result 

|2^{/9Vft?;} • H| <e(ll\/pVftn|||2(r^) +/3“|e“|||[Vftwllli2(e-) +H|e+|||[[Vft'f^llli2(e+)) 

i^‘^~^^)P nr' tiii2 11 [t tiii2 

^ ^ IIHIlL2(e-) H Xl+I IIHIlL2(g+). 


e e 
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Summing over all edges e G we get 

|2E fiP^hv}-Ml 


eeSl 




+ E 


(C^ + i)p- 


e\e 


[[^lllL 2 (e-) + 


(C 2 + l)pH 


e C 


^lllL2(e+)- 


Finally, the result follows by choosing e = 1/2 and then choosing 7 = - - + 1. 


□ 


4.2. Continuity of bilinear form. Next we prove continuity of the bilinear form. To do this 
we define the augmented norm 


(4.5) ||n||^ = ||n||^+E(^ U^hv} ■ n\\l2^^-^\e \ + P^\\{Vhv} ■ n\\l2(^^+^\e+\). 


Lemma 8. (Continuity) Suppose that v, w gV. Then, there exists a constant C > 0, indepen¬ 
dent of h, p~ and p'^, such that 

ah{w,v) < C lltcllvu Iblliu- 
Additionally, if w G V and v G Vh we have 


(4.6) 


ah{w,v) < C ||tc||vu ||u||\/- 


Proof. We give a sketch of the proof by bounding each term of ah{-, ■) in (2.6) separately. The 
first term can easily be bounded by Cauchy-Schwarz inequality 

/ pVhV-VhW < \\^"^hv\\L-^(^n)\\'/p^hw\\L 2 (^y 

Jn 

The first part of the second term can be written as 


f {pVhv} • M = [ {p Vftu} • M + [ • M- 

Je Je~ ^e+ 

Using the Cauchy-Schwarz inequality one has 

J^^{p^Vhv}lwj < • n||i2(e±)) -^^||HllL2(e±). 

Let e = Int(c?Ti n CIT 2 ). If n G V/ then Vhv\j,± is constant for each z = 1, 2. Then, we can use 
Lemma to get 

yi^\/p^ll{Vhu}-n||i2(e±) < ^A^\^/p^\\V^hvj\\L2{e±)+C^/p^{\\Vhv\\L2(^Tf) + \\'^hv\\L2(^Tf)) 
where C here depends on 6 . The third term can be bounded by 

[ < P'^IIMIlL2(e±)IINIU2(e±). 

J 

Finally, the fourth term can be bounded by Cauchy-Schwarz inequality 

[ < p^\\lVhvjh2^,±)\\lVhwjh2^,±y 

J 


The proof is complete summing over the edges, using finite overlapping of the elements associated 
to the edges and using arithmetic-geometric mean inequality. □ 
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5. A PRIORI ERROR ESTIMATES 
The purpose of this section is to prove a priori error estimates for the method defined in 


(2.7)-(2.6) in the energy and norm. 
5.1. Energy error estimates. 


Theorem 1. (A priori energy error estimate) Let u be the solution to problem (1.1) and Uh € Vh 
be its finite element approximation solution of (2.7). Then, there exists C > 0, independent of 
h, p~ and p'^, such that 

\\u-Uh\\v < Ch{l+n) (^^/p-{\\Du\\l 2 (q_-) + \\D‘^u\\l 2 (^q_-)) + \4^(llL2(n+) + \\D^u\\L2{n+))^ 
The constant C depends on Ce from (3.3). 

Proof. Note that ||m — Uh\\v < ||^ ~ Ihu\\v + \\IhU — Uh\\v- Since dh ■= IhU — Uh £ 14, coercivity 
of ahi-,-) (see ( |T^ ) gives 

(5.1) c\\dh\\v < ah{dh,dh) = ahihu - u,dh) + ah{u - Uh,dh). 


Using (4.6) we obtain, for e > 0 
(5.2) 


e„ , „2 1 


ahihu - u,dh) < C* ||4u -'u||iv||4||y < (7 |^-||4||y + —||4u - u||^ 

Choosing e sufficiently small we have 

(5.3) IM/illy ^ C \\IhU — uW"^ah{u — Uh,dh). 

The proof of the theorem follows by checking below Lemma[To| (approximation error in VU-norm) 
and Lemma 12 (inconsistency error). □ 

We next prove the missing parts of Theorem We need an estimate for interpolant 4 (see 


(3.5)) in the augmented norm W, and to do so, we first need to prove a trace inequality that 
goes from a part of boundary to the interior of the domain. 

Lemma 9. Let T G and let e be an edge of T. Suppose that w G U and 

le^l > 0. Then, there exists a constant such that 


1 


w 


Ii2(e±) < ^ (^fi2\Ml2^^±) + \\DM\l2(^^±)+hUD‘^w\\l2^T±)^ ■ 


Proof. Consider the case where enT 7 ^ 0. Let ei be an interior edge of ujt, contained in cofi and 
connected by one node to e~ . Edge ei exists thanks to the r-tubular neighborhood assumption 
with 2h < r. Note that cLt < \ci\ for a constant c that only depends on the shape regularity of 
the mesh. Using the fundamental theorem of calculus we can show (see Lemma 3 of [TO] for a 
similar two dimensional result) 

|g3|ll^lli2(e-) < C{\e I ||Du;||| 2 (e-) + ^||'W^|li 2 (ei) +/^T||T>'«;||| 2 (ei))- 

Noting that > Ch"^, and using standard trace inequalities gives the result. The exact same 
argument would apply for e^. If e n T = 0, e“ = e we apply trace inequality from e to uf. □ 

Lemma 10. (Best Approximation Error Estimate) It holds, 

\\IhU—u\\w < C h{l + k) {\\Du\\ e 2 (^^-) + ||-D^rt| 42 (Q-)) + \/^(||L>u| 42 (q+) + \\D‘^u 
H ere the constant C depends on the constant Ce from 
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Proof. We bound each term of 
easily have 


Ihu\\w: see (4.5) and (4.1), separately. Using Lemmawe 


||\/pV/i(tt < C (^p \\D u^\\^2f^rp-j + p'^\\D u~E\\i^2^rp'j^ 

T£Th\Tf 

Terr 

< C {1 + k)^ ^p {\\Du\\j^2(^Q--j + \\D‘^u\\j^2(^Q-)) + p’''(||I1r||^2(q+) + \\D‘^u\\j^2(^Q+-^)J ■ 


In the last inequality we used (3.3). We also have used finite overlapping of the sets lot- We 
next estimate the term 


ee£[ 


Suppose that e = Int(clTi n CIT 2 ). Then, using Lemma|^we have 


1^ II - hv\ Ili2(,-) +1^111^- 4 nl ||i2(,+) 

^ C* 1^(11“ “ -^Ti'*^lli2(e-) + 11'“' “ -^r2^^|li2(e-)) + C-^{\\u - Iti |||2(e+) + \\u - -^T2 Ili2(e+)) 


< 


E f - ^T^iii2(,-) + mu - )+4.ii^'(^ - ^r«)iii2(,-)) 

i=i \ r ‘ * » / 

■^'L'“IIl2(^^+j + ll-C'('“ - -^rJi)||^2(^+j + ^tJI-D^(?i - '^ri'w)||^2(^+j^ 


1 

— 

2=1 \ 


+ Cp'^ ^1 —lilt - 


If we now apply Lemma and use (3.3), we get 


E - 2i.“llli-(=-) + ^111“ - Ai*lili.,,+,) 


sg£F 


< Ch'^{l + nf [j) (||ZlR||^2(f2-) + Il-C^^'“lli2(f2-)) + P^(II-^^IIl2(q+) + ||L>^tt||^2(f2+))) 


We next bound 


E 1 ^ I IIP(u- 4 u)l||^ 2 (g-) +p+|e+| ||p(R- 4 R)l||^ 2 (e+)J . 
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Using a trace inequality to obtain 

|2 


P |e I ||[[Vh(n-4n)l||^2(g_)+p+|e+| ||[V/*(n-4n)l||^2(g+) 

< p-\e\ ||[V,(n - 4n)l||i2(,-) + p+|e| ||[V;,(n - 4«)llli2(e+) 

< P“|e| {\\Diu - /Tili)|li2(e-) + \\D{u - /r2^^)|li2(e-)) 

+ P'''kl “ -^Ti^)lli2(e+) + \\Diu — 42^i)|li2(g+)^ 


2=1 






- >, + 11 (ti — u 




+ Y. ) + h^TMD\u - ) 

i=l ^ 

Again, if we apply Lemma and (3.3) we obtain 

E (p“|e"llllV;,4-4«)llli2(g-) +p+|e+|||[V,4-4n)l||i2(g+) 


< C/i^(l + k )^ (||L)tt||^2(f^-) + ||L)^tt||^2(f^-)) + p’''(||L>rt||^2(Q+) + ||L)^tt||^2(f^+))^ 

Finally, the last term is estimated by 

Y i\^~\P~\\{^h{u - Ihu)} ■ n\\l 2 ^^-) + |e+|p+||{Vh(n - Ihu)} ■ n||^ 2 (g+)) 

T&Sh 

< Ch?{l + [p (||D'u||^2(f^-) + ||L>^'a||^2(f^-)) + p’^(||L>'u||^2(q+) + ||L>^n||^2(f^+))^ 
exactly in the same way as it was done above. 

In order to bound the inconsistency term we need a technical lemma. 


□ 


Lemma 11. Let T G and enumerate the three edges of T by 61 , 62 , 63 . Then, there exists a 
constant m such that 

iTrl < m max \e~\ and |Tr| < m max |6i*~|. 
i=l,2,3 j=l,2,3 

The constant m is independent of r (see Lemma^, depends only on the shape regularity of the 
triangulation. 


Proof. See Appendix |A.2[ 


□ 


Now we are able to establish the inconsistency error estimate. The constant will be indepen¬ 
dent of the contrast of the coefficients p~ and p'^. 


Lemma 12. (Inconsistency Error Estimate) Let u be the solution to (1.1) and Uh € Vh be its 

finite element approximation solution to (2.7). Then for any dh G 14, it holds 

(5.4) 

ah{u - Uh,dh) < Chn(^^ffr {\\Du\\l 2 (^q-) + /i||L)^m| 42 (q-)) + ^/fFh\\D'^u\\L 2 (^Q+)'j \\dh\\v- 
The constant C depends on the constant Ce from 
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Proof. First note that ah{u — Uh, dh) = ah{u, dh) — (/, dh). The inconsistency term is then given 
by 

ah{u-Uh,dh) =yZ iP~Dn-u~)[dh], 

JTr 


Terr 


where we have used that p~Dn-u~ = p'^D^+u^ on F. We also have used the notation \dh] = 
dj^ — d~f. Let Lt be the tangent line to F at xq and note that [d^] is linear function on Lt 
that vanishes xq as well as its first derivative. Hence, [dh] = 0 on the line Lt- Since the 
distance between this line and the curve Tr is O(s^), where we set s = |Tr|, we have that 
|[(i/i](x)| < C11X" 11 2,00 for every x G Tp. Hence, we have 

/ {P~Dn-U-)[dh] < C ^^s‘^\\^/p^D^-Uf;\\L2(^Tr)^/^\\D[dh]\\L^Tr)■ 

■JT^ 

For the moment, let us assume that following inequality holds 
(5.5) ^/p^^/s\\D[dh]\\L^T^) < CMt, 

where 

Mt = '^ (^/\e^\\^/f^lVhdhj\\L•2(e+) + —7F^\\^/^dh\\L2(e+) + \\^/^'^hdh\\L2(^^+)j 

+ X] ( V\^~\\\y^Whdh}\\L'^{e-) H-77=^ll'\/^'^/illL2(e-) + I • 

eCdT \ V r I / 


Then, 


[ ip D^-U )[dh] < C KS^/^\\^/p^DuJ^\\L 2 ^T^)MT. 

J 7p 


Letting Bg be the ball of radius s centered at xq we can use the trace inequality to get 

y/s 


Hence, 


/Tr 


[p D^-u )[4] < Csdlv^ Duj^\\l 2 (^Bs) + 4Vf^ 11 2,2(5,)) Mr. 


We see that (5.4) follows after using that 

< C\\dh\\l, 

Terr 


the inequality (3.3) and the fact that s < C Ht (see Lemma A.l). 


In order to complete the proof we need to prove ( |5.5[ ). Firstly, by using the triangle inequality 
we have 

\/^\/s||L)[4]||T2(rp) < \/s (||L)4 ||2,2(rr) + ||T)(i)J)||2,2(Tr)) ■ 

Then by Lemma 11, there exists an edge e of T such that 

(5.6) s < m|e'''|. 

Now let e = Int(cIT n dK), fox a K G Th- Using that Dd'^ is constant we get 


l|T>'^/l(llL2(rr) — \fs\Dd'l\ — 


f\ X 


^l|T»d)|;||L2(e+) < Vm\\Ddl\\L2< 


e+)- 
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According to Lemma 

(5.7) |e+|2 < 0max{|r+|,|A:+|}. 

First, suppose that |r+| = max{|T+|, |A+|}. Then, 

l|Vfcd+||i2(e+) = -^Q||D(i+||i,2(r+). 

Hence, 

l|V/.<k 2 (,+) < ^\\Dd+h 2 ^T+y 

V'5 

On the other hand, if |Ar'’'| = max{|T^|, |iir'''|} then we have 

l|V/i^^/l(llL2(e+) < \\D{dl\T - d'l\K)\\L2{e+) + -^-^\\Ddl\\L2(^K+y 

We write D{dh\T — dhlx) as its normal part and tangential part, then, use an inverse inequality 
for the tangential part on e"*“, and have 

||V/,(d;J)|r - d/J)k)||L2(e+) < ||[Vft4]]||i2(e+) + |^||[4]]||L2(e+), 

therefore, 

\/f)m O 

\\'^hd'l\\L'2(e+) < l|[[V/j4]]||x,2(e+) + —^ ||x, 2 (^+) + |^ || [41 ||L 2 (e+). 

Combining the above inequalities we obtain 


•/F~40d+U.^Tr) < y7^\Tf|l|Vi<iJ||i«„+, + 

+ Cm\/7(||vTr>AllL2(r+) + llvTr’AllL2(7^+)). 

where we have used p~ < p'^, so 

\/ p-s\\Dd'l\\L2pp-^-^ - ^ X/ ('\/"^|e+|\/4ll[^/j41||L2(g+) + l|[41||L2(g+)) 

eddT vl® I 

+ m\/0||v^Vfe4||^2(^+). 

Using the same argument we can prove 


{p s)^/2 ||L»4 42 (rj,) <CY^ (\/m|e |v^ ||[L>41l|L2(e-) + '^r. -|^ ll[41l|L2(e-)) 

ecar Vle I 

+ Cmv/^11 Vw Vh411 i2 . 

Combining the two last inequalities we obtain ( |5.5[ ). □ 

Now we would like to state a corollary of Theorem We first need some regularity results. 
We start with a standard energy estimate. 


Proposition 1. Let u solve (1.1), then the following energy estimate holds 


(5.8) 





II/IIl 2 (o) 


where C is independent of p^. 
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At this point we recall that we are assuming that p < /o^. In fact, a better energy estimate 
holds as we prove in the next proposition. 


Proposition 2. Let u solve (1.1), then the following energy estimate holds 

(5.9) p'^\\Du\\l2(q+) + p~\\Du\\L2i^fi-^ < C ||/||i2(f^). 

The constant C could depend on the geometry of the subdomains 

Proof. We prove this estimate in the case that is the inclusion (i.e. Sn"*" does not intersect 
dLl). The proof of the other case is similar. First, from the previous proposition we have 

(5.10) p \\Du\\l2^^-'^ < (^11/111,2(^2). 

Let w = u~^dx. Therefore, using Poincare’s inequality we have 

(5.11) ll^llL2(n+) < C\\Du\\L2(^Q+y 

By means of an extension result (see Lemma 6.37 in |lljL there exists v G Hq{LI) such that 
u|q+ = m|Q+ such that 

(5.12) — ^ll'“^llHi(r2+)- 

Applying the variational formulation of problem (O) we have that 

/ p~Du ■ Dv + / p'^Du ■ Dvdx = / fv. 

Jo.— Jft 


Multiplying by p'^ it follows that 

r,,.ll 2 


in+ 


p'^Du ■ Dvdx = p~^ fv — p'^ 


p Du ■ Dv 


IQ- 


IIp ^^IIl2(c+) = p 

Therefore, we have 

||/9+T>u|||2(f2+) < \\f\\D{n)\\p'^v\\D{Q) + ll/5".DR||i2(n-)||/9+L>u||i2(f2). 

By ( |5.12 ) we have 

(5.13) ||p+Z)R||^2(f2+) < (^(II/IIl2(0) + Up .^^llL2(0-))P^||'W^||iLi(0+). 

Hence, using ( 5.11[ ) we get 

(5.14) ||p+L>n||^2 (Q+) < (1^(||/||l2(q) + Up Du\\L 2 ^Q-^)\\p'^Du\\L 2 (^Q+y 

The proof is complete after applying ( |5.10[ ). 

We will also need to state and regularity estimates which can essentially be proved if one 
combines results in [8] and m- Here we give an alternative argument using the results in [T^ 
and classical regularity theory. 

Proposition 3. In addition to the assumptions already made in this article, assume further 
that Ll is convex. Let u solve ([ 13 . Then, the following regularity estimate holds 

(5.15) p'^\\D'^u\\i2(^q+) + p \\D‘^u\\i2(^q -) ^ Creg II/IIl2(o), 
where C is independent of p^. 

The constant C^eg could depend on the geometry of the subdomains 


□ 
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Proof. We consider two cases: 

Case 1: is the inclusion (i.e. dQ~ H 917 = 0) 

In this case, (5.15) follows from (2.35) in [T^ . 

Case 2: n"*" is the inclusion 
From (2.36) in [19] we have 

(5-16) 1^11^-2(11-) < C” ll/llL2(n)- 

Let u'^ = u~^dx. Let w = u~^ — u~^. Note that w satisfies 

—= / in 17“'“, 

p'^Vw ■ n = p~^Vu~^ ■ n on 917“'“. 

From a standard regularity result (see (2.3.3.1) in [T3j) we have 

P'^\\w\\min+) < C'(II/IIl2(q+) ||p“^Vu“^ • n\\Hi/2(^QQ+-^ -h /3 +||r^|Ui(o+)). 
Using the jump condition we have 

||p^Vn“'“ • = Up Vu ■ n\\fji/ 2 (^QQ+y 

By a standard trace inequality we have 

||p"Vu" • n\\fji/2(^Q^+) < p"||w“||/i2(o-)- 


Which by (5.16) gives 

||p“^Vn+ • n\\jji/ 2 (^g^+) < C Wfh^n)- 

The estimate 

- ^ll/llR2(fi) 

follows from Poincare’s inequality and Proposition!^ Hence, we have 

P^\\w\\H^(n+) < C* ll/llL2(n)- 

The result now follows after noting that \\D‘^u\\l^[q+) < ||r’||/i 2 (q+). 

Combining Theorem and the previous propositions we have the following corollary. 


□ 


Corollary 1. Let u be the solution to (1.1) and Uh G 14 be its finite element approximation. If 
17 is convex, we have 

Ch{l + n). 


(5.17) 


\u-Uh\\v < 


^—ll/llL2(n)- 


The constant C depends on the constant Ce from (3.3) and the constant Creg from (5.15). 
5.2. L^-Error Estimates. We now prove an error estimate using a duality argument. 


Theorem 2. Let u be the solution to problem (1.1) and Uh ^ be its finite element approxima¬ 
tion solution to (2.7). Assume that 17 is convex. Then, there exists a constant C > 0 independent 
of h, p~ and p+, such that 

C /i^(l + k)^ 

(5-18) \\u-Uh\\L^(n) < --II/IIl2(o)- 

The constant C depends on Ce from (|3.3[) and Creg from (|5.15[). 
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Proof. Let (f be 

a solution of the problem 


(5.19a) 

-p^Acj)^ = {u- Uh)^ 

in 

(5.19b) 

-c- 

II 

o 

on 5H 

(5.19c) 

o 

II 

on T, 

(5.19d) 

[pDn(f>] = 0 

on T. 

We have 




ll« -'“/illizm) = / {u-Uh)u- / {u - Uh)uh = ah{(p,u) - ah{4>h,Uh), 

JQ Jn 

where (l)h is the finite element approximation of (j). Therefore, we see that 
(5.20) = ahi(i),u-Uh) + tthicj) - (ph,Uh) 

^ ^h) T 4^h) T 4‘ht ^/i) 

^ ^/i) T ^h4^) T '^ht 

T 4^h} T ^h^)‘ 

Using continuity of the bilinear form we have 

ahi4> - 4>h,u - Uh) < C\\4i - (phWw \\u - UhWw- 

Using the triangle inequality we have 

11?^ - "W/ilIvU < lltt --f/iR||vu + llAw - Jthllvu- 
It is not difficult to show that 


\\IhU-Uh\\w < C\\IhU-Uh\\v- 


Hence, 


\\u — Uh\\w < \\u — Ihu\\w + \\u — Uh\\v- 

Now using Theorem 10, Theorem ( |5.15 ) and (5.8) we get 

Ch{l + K) 


\u — UhWw < 


■||/llL 2 (n)- 


ill ^ h(l -)- /^).. .. 

- (ph\\w < - -;= — \\u-Uh\\L^(n)^ 


Similarly, 


and hence we have a bound for the first term in ( 5.20| ) 

(J }{^(\ -\- ^^)2 

ah{<t> - (t^h.U - Uh) < -z- \\f\\L 2 {Q)\\u-Uh\\L^(n). 


Using Lemma 12, (5.15) and (5.8) we have 


ah{u- UhAh- Ih4>) < ^ II/IIl2(o)II4</> - (t>h\\v, 


which implies the following bound for the second term in (5.20) 

ah{u - UhAh - Ih(l>) < ||/||^2(i^)||u - Ufe||^2(i^). 
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Analogously, for the fourth term in (5.20)we have 


(j 

0^hi4’ ~ 4^hjUh — IhU) < ^ ll/llL2(f2) 11'^ ~llL2(f2). 


For the third term in (|5.20|) we have 
(5.21) 


ah{u- Uh,Ih4') = ip D^-u )[4(/>]. 

Using the Cauchy-Schwarz inequality we obtain 


1/2 


ahiu-Uh,Ih(t>) < \\\fp~ Dn-U ||L2(r)\4r ^ WihMh^Tr) 


For the first term in (5.2) we apply a trace inequality to obtain 


C 


\\y/p~ DnU ||L2(r) < C^ffT (||F>^u||l2 (0-) + \\Du\\l2(^q-)) < \\f\\L2^Qy 


where we used (5.15), (5.8) and the fact that p~ < p'^. Now for the second term in (5.2), we 
note that = 0 on Lt and that Lt is at most distance 0{s^) from Tr where s = |Tr|. 

Therefore we can use, Taylor’s theorem to show that 

[Ih(t>]{x) < C\\X''\\L^s^\D[hcl>]\{x) VxGTr. 

Hence, 

l|[-^h9^]||L2(rr) — C^x'^\\D[Ih4>]\\L‘2(TY) — ^'^hT\\^[^h4>]\\L^{Tr)^ 
where we used that s < C hx- Consequently, adding and subtracting D[<p] 

1/2 / 1/2 

\t&tF 


< 


Ck Y1 hUD[h(^-mhiTr) E hUmWUTr) 




,Terr 


Observe that 

h‘^\\^[4>]\\L^{T) < Ch^{\\D(t)^\\x2(Y-) + WDcf) ||L 2 (r)). 

Application of trace inequality gives 
\ 1/2 

E ll-^HIli2(Tr) I - (II'^^<?^IIl2(q-) + \\D4>\\l2{p,-) + \\D‘^<t>\\L'^{n+) + \\D(j)\\L2(fi+)) 

For the other term we use that Ih(p\T = It^' S'Hd use a trace inequality to bound 

\\D[Ih(p — 4’]\\L‘^(Tr) = \\D[It 4> — 4’]\\L^{Tr) 


From (3.6) and (3.7) (using that p < p'^) we obtain 

l| 0 (W - + l| 0 (W - ) < cfcT(l|oVlli.,„-, + ,). 
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Therefore, 


1/2 


hUD[h^-mhiTr) 1 < C/iVPV||l2(C-) + P</>IIl2(0-)) 

+ C'/i^(||Z/^(/)||^2(Q+) + \\D4'\\L2{n+))- 




Hence, using the regularity result (5.15) and (5.8) we have 




< — =\\u - Uh\\L2(n)- 


\4^ 


Thus, we obtain the bound for the third term in (5.20) 

Ch?K 

ahiu- Uh,Ih(p) < -^^\\f\\L^n)\\u-Uh\\L2(n)- 

In a similar fashion we can prove 

Ch?K, 

ahi(i>- 4>hJhu) < -^II/IIl2(o)II^^ -''^h||L2(o)- 

The proof is complete after combining the above inequalities for the terms in (5.20). 

6. Extensions and final remarks 


□ 


6.1. Extension to three dimensions. We now show that the space S^{T), introduced in 
(2.4), can easily be extended in three-dimensions. We consider the problem 0 where H is 
a three-dimensional domain and T is a simple, closed, surface. Let now 7h he a simplicial 
triangulation of H. Let be all the faces of the mesh Th- Finally, let be the set of all faces 
that belong to a tetrahedron that intersects T. 

We now define the local hnite element space. Let T £ Th he a tetrahedron. Let xq be a hxed 
point on T n T. Let tIq be the outward pointing unit vector normal to at xq. Let be 

such that forms an orthonormal system. 

Given v G P^(r+) there exists a unique T(u) G P^(r“) satisfying 

T{v){xo) = v{xo), 

(■^t+^(^)) (^o) = (®o), 

(^s+^(^)) (®o) = (Ds+v) (xo), 

P~ (^n+^(^)) (®o) = /5+ {Dr^+v) (xo)- 

Given T G T^ and for each v G P^(T^) we can consider the unique corresponding function 


Giv) = 


Tiv) 


in r+, 
in T~. 


Let span {vi,V 2 ,V 3 ,Vi} be a basis for P^(T) restricted to T+. Then we define the local finite 
element space 

= |span{G(ui),G(u2),G(u3),G(u4)} , if T G 7^^ 


pi(r) 


, if^GrAr^ 
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Consequently, the global finite element space is given by 

Vh '■= : v\t G 5^(r), VT G Thi V is continuous across all faces in . 

We believe that can be extend and analyzed to the three-dimensional case, where now 
e± = e n 0=^ are pieces of a face e of a tetrahedra, and considering appropriate stabilization 
parameters. This could be a subject for a future work. 


6.2. Alternative Local Spaces. An alternative approach to enforce weak continuity of our 
local space S^{T) is normally used in the literature (see [S]). Instead of enforcing continuity 
of function and tangential derivative at xq one imposes continuity at two distinct points xi,X 2 - 
More precisely, define xi,X 2 be the two points in T that intersect dT. Then, one can define T(u) 
(in contrast to the definition in Lemma by 

T{v){xi) = v{xi) for i = 1,2 
P~iD^+T{v)){xo) = p+{D^+v){xo). 

Subsequently, Itu is now defined by 

{l:^u){xi) := (JrUg)(xi) =: (7^u)(xi) for i = 1,2 
p-{D^+I^u){xo) := p+{D^+Jtu^){xo) =: p-^{D^+l:};u){xo). 

Defining the finite element method using these local spaces, we can prove all the a-priori esti¬ 
mates above. 

Another alternative is to enforce the matching conditions by averaging on Tp. More precisely, 
one can redefine T(u) in Lemma by 



or alternatively, one can replace the second equation by 



Df+T{v) ds 



D^+ vds. 


Even though their numerical implementation is more complicated and require numerical inte¬ 
grations, these alternatives have the advantage that the analysis is shorter especially for the 
consistency error, and also can be formulated using Lagrange multipliers. 


6.3. Extension to Cartesian grids. We note that the analysis and methods developed here 
can be easily extended to Cartesian grids and elements. One possibility would be to use 
same local spaces S^{T) defined before for quadrangular elements T G that is, one has 
three degrees of freedom for T G and four degrees otherwise. In case one wants to have four 
degrees of freedom for T G , let the space S^{T) then be a piecewise bilinear function under 
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the coordinates no, to and impose 

T(r;)(xo) 

{D^+T{v))ixo) 

P~iDr,+ '^{v)){xo) 


n(xo) 

{D^+v){xo) 

P'^iDr,+ v)ixo) 


6.4. Alternative Bilinear forms. In this section we discuss alternative bilinear forms. We 
will point out what are the theoretical difficulties in analyzing these alternative methods. When 
we provide numerical experiments in the following section we will also see that, although we 
cannot prove stability for some methods, they sometimes do well experimentally in some of the 
norms. 

Firstly, formulation without flux stabilization has been used for example in |2T]. Methods 


( |6.1[ ) and (6.2) experiment this direction. 

(6.1) ah{w,v) := J pVhW ■ VhV - j^{{pVhv} ■ {wj + {pVhw} ■ {vj) 




+ 7 X] I / P HJI • + 




p iiw^ii • m 


( 6 . 2 ) 


ah{w,v) := / pVhW-VhV-'^ / ({pV^u} • [u;]] + {pV/jU;} • H) 






p + 


P • FJ 


The problem with these methods is that we cannot establish neither the optimal inconsistency 
error nor the coercivity independent of the contrast and the mesh. 

We note that methods (6.1) and (6.2) seem to do well numerically in the L? and energy norms, 
however, they are mesh dependent in the L°° norm and they do not converge in the weighted 
norm. 

The method that seems to do the best numerically, slightly better than the proposed method, 
is the following one: 


(6.3) 


ah{w,v) := / pVhW-VhV-'^ / ({pV/^u} • [u;]] + {pV/iw} • [uj) 


e&el 




p”H • H + 


p ' • lluj 


IF ^ |e| p [V/iu]] [Vfeu;]] + p'^l'^hvj [V/iU;]] 




The difference between this method and the one analyzed in this paper is that we are using 
a stronger flux stabilization (i.e. we replace |e^| by |e| and we introduce a flux stabilization 
parameter jf)- It is not difficult to see, that we can prove all the error estimates contained in 
this paper for this method. In particular, the coercivity of the bilinear form is obvious since we 
are adding even more stabilization. Clearly, now we would have to redefine our V and W norms 
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but the approximation properties will still hold. In summary, Theorem and Corollaries and 
[2] hold for this formulation. 

A very natural question would be if we can penalize the jump terms with l/|e| instead of 
l/|e^| and get a stable and optimally convergent method independent of contrast. The answer 
is yes, as long as we penalize the jump of the full gradient instead of the flux. 

(6.4) ahiw, v) := pVhW ■ {{p'^hv} ■ M + {p^hw} ■ H) 




1 




+ ^ E y 1 / P~M-M + I P'^M • 


7f |e| p [Vu]^ • [Vin]^ + [Vu]^ • [Vw]^^ 

where the jumps in the last line are defined by 


[Vn]^ = Vu (g) n +Vu’' 


i n 


and the symbol ® denotes the outer product between vectors. We note that in the coercivity 
and inconsistency error analysis proofs, we have used the inverse inequality \e~\ \\Dnv\\L‘2(^e-) — 
C II (where T~ is a triangle with edge e, and similar result for e^) in particular in the 

estimation of the left-hand side of (4.4). For method (6.4) we want to avoid using this inverse 
estimate and the factor l/|e“| (also l/|e'’'|). Let x be the vertex of e belonging to 0“. Using the 
ideas in the proof of Lemmaj^ we can find a sequence of elements {T2, Ts, • • • , T/v} in the patch 
of X, and common edges {e2, 63, • • • , eM-i} with e* = T, n Tj+i and e~ := et n 0“, such that, 
|e“| < (T|e2 I < o'^lejj’l < = (T^~^\e\ where a a positive constant bounded uniformly 

from below by zero and N is bounded depending on the shape regularity of the mesh. Then, we 


bound the left-hand side of (4.4) by using recursively 


and 


and then use that 


Hence, 


||Vu|ri_i ||^2(e-) ^ 2 ^11 [Vu]g, 11^2(6-) + ll^^lTllL2(e-)) ) 
l|Vi;|rJ|^2(,p < l/cT||Vn|Tj|^2(,-^^), 

|e| l|Vu||i2(,-) < c |e| II [Vu]^ + - \\Vv\\Y^.y 


i=2 

6.5. Final remarks. In this article rigorous error estimates independent of contrast were de¬ 
rived. However, many interesting research questions remain. Firstly, we kept track, as much 
as possible, of how the constants depend on the geometry (e.g. curvature). However, there 
are two constants, Ce and Creg, that we did not investigate in detail how they scale with geo¬ 
metric quantities such as maximum curvature and the radius r of the tubular neighborhood. 
We believe it is possible to prove a bound for Ce in terms of geometric quantities, but for 
the sake of simplicity we did not investigate it here. In addition, the regularity constant Creg 


appearing in 5.15 depends on the geometry. Addressing these issues will lead to error estimates 
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that are completely explicit on their dependence on the curvature and radius r of the tubular 
neighborhood. 

Secondly, our method and error estimates do not consider problems with high curvature. 
An interesting line of research would be to investigate problems with interfaces of arbitrary 
curvature. In this direction, perhaps a combination of the method presented here and the 
multicale method by Chu et al. in [S] could be a possible approach. 

Thirdly, rigorous error analysis for higher order Immersed Finite Element methods (fe > 1) 
remains open. It would be interesting to study this, in particular, for high-contrast problems. 

Finally, we would like to study the sharpness of our Assumptions (l)-(3). For instance, we 
observe that our method is still well-defined without Assumption (3). 


7. Numerical examples 

In this section we explore the properties of the methods presented in sections above applied to 


the two dimensional interface problem (1.1). In particular, we are interested in the computation 


of the following errors and their respective estimated order of convergence 



'■= \W - Uh\\L'^{Q.)^ 

ef 



■= 11 ,, (u-«/,)! 1^2 (Q), 

1,00 

■= \\VP^h{u-Uh)\\L^(n), 

4 

:= \\pSIh{u-Uh)\\L^(n), 

- 1,00 

:= ||/9V/i(u - Uh)llioo(f^), 

Tl.OO 

h 

:= \\pDn{u-Uh)\\L^(T), 

- 1,00 

:= \\p^h{u — 'W/i) L°°(r2\(U{r:Te4'"})) 


e.o.c. := 


^ogiehi.Jeh 


\og{hi+i/hi) ■ 

Computation (or approximation) of the L°° norms are performed evaluating the error at the 
following set of points: if an element does not intersect the interface then the set of points 
are the nodes and the centroid of the element, if the element does intersect the interface we 
evaluate at the nodes and at the two points intersections of the interface with the edges of the 
element. We observe that the errors corresponding to the norm and weighted with 
semi-norm correspond to the only results proved in this paper. Theorem and respectively. 
We expect optimal convergence, second order in the norm and first order in the weighted 
semi-norm. We also compute the analogous in the norm and W^’°° weighted with 
semi-norm. In addition, we compute the errors in the IF^’^ and semi-norm weighted with 

p. Note that the estimate for the interpolation error is also achieved in this semi-norm. Indeed, 
this is a consequence of Lemma and Proposition and 

||pVfe(u-4u)||i2(t5) < C\\f\\L2(n). 

A revealing result of our experiments is that the ratio of convergence of the error is optimal 
when we compute using the triangles non intersecting the interface. Error illustrates this 
observation. Finally, error is a standard error for interface problems, illustrating the 

approximation of the normal derivative on the interface. 


The experiment presented below shows that method (2.7)-(6.3) produces the best results. 


This method (and all the others) does not seem optimal for the error weighted with p. 

However, it is optimal when we do not consider the elements in . We highlight that this 
result was not remotely addressed in this paper, and this kind of estimates appear to be more 


difficult. For the method that we analyzed (2.7)-(2.6) we observe an uncertain behavior in the 
error 6^’°° for the last mesh, not achieving the optimal convergence as in the previous method. 


We also present tables for one of the method without flux stabilization, method (2.7)-(6.2). We 
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observe non convergence so far for the error e^’°° and in the last mesh and a deterioration in 
the norm. We also do not have the optimal convergence for the normal flux error 

In our numerical experiment we consider the two dimensional domain = (—1,1)^ with the 
immersed interface T = {x ^ Q ■. x\ + = E?}. We define 11 “ := {x ^ Q. : x\ + x\ < and 

= Q\{Q~ U r). Example |(l)|considers the following exact solution 


u{x) 




+ - 
^ 3“ 1 p- 


if ® G Q , 
if a; G 0 +, 


where R = y/x^R X 2 and a = 2. For example ( 1 ) we have E = 917“. Similar results were 
obtained for the case E = Sn"*". We provide plots of the approximate solution for both cases. 

Finite element uniform triangular meshes Th non matching the interface were used. In the 
tables we compute with h = for I = 1, ..., 7. 

All the computations were performed in MATLAB, including solution of the linear system by 
means of “\”. 

( 1 ) Case E = 917 , p'*' = 10^ and p =1. Tables 0 and show the results obtained by 
method (2.7)-(6.3) with stabilization parameters 7 = 10 and 'jp = 10. Note that this 
method has a stronger flux stabilization that the method analyzed in the paper. 


1 

pO 

e.o.c. 

pOO 

e.o.c. 


e.o.c. 

1,00 

e.o.c. 

1 

8.2e-3 


2.5e-2 


l.le-1 


3.7e-l 


2 

1.7e-3 

2.28 

5.7e-3 

2.12 

4.4e-2 

1.30 

2.1e-l 

0.86 

3 

2.7e-4 

2.63 

1.3e-3 

2.19 

1.8e-2 

1.29 

9.7e-2 

1.08 

4 

4.6e-5 

2.57 

3.2e-4 

1.97 

8.3e-3 

1.12 

5.2e-2 

0.90 

5 

9.0e-6 

2.34 

7.2e-5 

2.15 

3.9e-3 

1.07 

2.5e-2 

1.08 

6 

2.0e-6 

2.19 

1.8e-5 

2.01 

1.9e-3 

1.03 

1.3e-2 

0.92 

7 

4.7e-7 

2.08 

4.6e-6 

1.94 

9.5e-4 

1.02 

6.8e-3 

0.94 


1 

4 

e.o.c. 

-1,00 

eE 

e.o.c. 

~l,oo 

e.o.c. 

n,oo 

e.o.c. 

1 

3.9e-l 


7.0e-l 


7.0e-l 


3.7e-l 


2 

1.6e-l 

1.32 

6.1e-l 

0.19 

4.0e-l 

0.81 

2.1e-l 

0.86 

3 

6.4e-2 

1.29 

1.9e-l 

1.68 

1.9e-l 

1.07 

9.7e-2 

1.08 

4 

2.9e-2 

1.15 

2.2e-l 

-0.20 

l.Oe-1 

0.91 

4.9e-2 

1.00 

5 

1.4e-2 

1.09 

1.4e-l 

0.65 

5.0e-2 

1.01 

2.5e-2 

0.98 

6 

6.6e-3 

1.04 

6.6e-2 

1.09 

2.5e-2 

0.99 

1.2e-2 

1.00 

7 

3.2e-3 

1.02 

5.1e-l 

-2.95 

1.3e-2 

0.98 

6.2e-3 

1.00 


Table 1. Example (1) 
10“^, using method (|2.7|) 


err 

7D-(|6.3D 


errors and convergence orders with 
with stabilization parameters 7 = 


p =1 and p+ = 
= 10 and 7 i? = 10 . 


The results in Table show optimal convergence for the errors e\ and validating 
the theoretical results Theorem and In addition we observe optimal convergence 
for the error e“ and e^’“. The error e\, weighted with p instead of ^/p, converges opti¬ 
mally. However, the rate of convergence of the error is not optimal. An interesting 
observation is that the error converges optimally, indicating that is only a couple of 
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-1,00 

e-h 

IQi 

2.3e-6 

6.5e-3 

2.8e-3 

10^ 

2.0e-6 

6.6e-3 

2.0e-3 

10^ 

2.0e-6 

6.6e-3 

1.9e-3 

10'‘ 

2.0e-6 

6.6e-3 

1.9e-3 

lO'^ 

2.0e-6 

6.6e-3 

1.9e-3 

10® 

2.0e-6 

6.6e-3 

1.9e-3 


Table 2. 


(2.7)-(6.3) 


Example (1) errors with p = 


with stabilization parameters 7 = 


1 and h = 2 ( 6 + 3 / 2 )^ 
10 and 7 ir = 10 . 


using method 


elements where the error does not converge. Another appealing feature of the method 
is the optimal order of convergence for the error . 

Table shows that the errors are independent of the contrast. We can observe that 
although we increase p'^ the errors remain constant, showing the contrast independency 
of our estimates. 

We present in Table the errors and convergence orders for the method analyzed in 
the paper (2.7p-(|2.6|). 


1 

pO 

e.o.c. 

„oo 

e.o.c. 


e.o.c. 

1,00 

e.o.c. 

1 

8.6e-3 


2.6e-2 


l.le-1 


3.7e-l 


2 

1.7e-3 

2.31 

6.0e-3 

2.14 

4.5e-2 

1.30 

2.0e-l 

0.90 

3 

2.8e-4 

2.63 

1.3e-3 

2.18 

1.8e-2 

1.31 

9.3e-2 

1.09 

4 

4.7e-5 

2.57 

3.4e-4 

1.97 

8.4e-3 

1.12 

5.5e-2 

0.76 

5 

9.2e-6 

2.35 

7.6e-5 

2.16 

4.0e-3 

1.08 

2.6e-2 

1.09 

6 

2.0e-6 

2.20 

1.9e-5 

2.01 

1.9e-3 

1.03 

1.4e-2 

0.91 

7 

4.7e-7 

2.09 

5.1e-6 

1.88 

9.5e-4 

1.02 

7.1e-2 

-2.37 


1 

-1 


-1,00 


-1,00 


71,00 



e.o.c. 


e.o.c. 


e.o.c. 


e.o.c. 

1 

4.0e-l 


7.4e-l 


7.4e-l 


3.7e-l 


2 

1.6e-l 

1.33 

2.5e-L0 

-1.79 

4.0e-l 

0.90 

2.0e-l 

0.90 

3 

6.4e-2 

1.31 

2.8e-l 

3.20 

1.9e-l 

1.07 

9.2e-2 

1.11 

4 

2.9e-2 

1.15 

1.6e-L0 

-2.56 

9.7e-2 

0.96 

4.5e-2 

1.01 

5 

1.4e-2 

1.09 

4.8e-l 

1.76 

4.7e-2 

1.04 

2.3e-2 

0.96 

6 

6.6e-3 

1.05 

4.8e-l 

0.01 

2.4e-2 

1.00 

l.le-2 

1.05 

7 

3.2e-3 

1.02 

7.1e-L0 

-3.88 

1.8e-2 

0.42 

5.9e-3 

0.93 


errors and convergence orders with p =1 and p^ = 


Table 3. Example (1) 

10“^, using method (2.7)-(2.6) with stabilization parameters 7 = 10. 


The results in Table show optimal convergence for the errors ej^ and validating 
the theoretical results Theorem and In addition we observe optimal convergence 
for the error e“. The error seems to converge optimal up to mesh 1 = 6, and then 
for the last mesh the rate can possibly be affected by the choice of the flux stabilization 
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parameter. As in the previous test the error ej^ converges optimally, however for the rest 
of the errors the convergence is not as clear as in the previous method. We suspect that 
this phenomena is related to the weights |e^| in the flux stabilization. 

Finally, Table displays error and convergence orders for the method without flux 
stabilization ( |2.7[ )-( [6^ . 


1 

pO 

e.o.c. 

„oo 

e.o.c. 

el 

e.o.c. 

1,00 

e.o.c. 

1 

1.7e-3 


7.0e-3 


5.8e-2 


2.2e-l 


2 

4.5e-3 

1.88 

2.9e-3 

1.29 

3.1e-2 

0.89 

1.6e-l 

0.51 

3 

3.6e-4 

0.31 

6.1e-3 

-1.08 

1.2e-l 

-1.89 

1.3e-Fl 

-6.37 

4 

2.8e-5 

3.69 

2.1e-4 

4.87 

7.7e-3 

3.92 

2.0e-l 

5.98 

5 

7.3e-6 

1.96 

5.0e-5 

2.05 

3.8e-3 

1.02 

1.5e-l 

0.40 

6 

1.7e-6 

2.06 

1.2e-5 

2.05 

1.9e-3 

1.01 

2.5e-2 

-0.70 

7 

4.5e-7 

1.95 

4.0e-6 

1.60 

9.5e-4 

0.99 

l.le-1 

0.23 


/ 

4 

e.o.c. 

—1,00 

e.o.c. 

~l,oo 

e.o.c. 

71,00 

e.o.c. 

1 

2.1e-l 


5.8e-l 


2.2e-l 


1.5e-l 


2 

1.3e-l 

0.69 

1.4e-Fl 

-4.58 

1.5e-l 

0.49 

l.le-1 

0.39 

3 

7.3e-l 

-2.46 

1.3e-F2 

0.10 

8.9e-l 

-2.5 

1.3e-Fl 

-6.83 

4 

2.8e-2 

4.70 

7.2e-F0 

0.84 

5.3e-2 

4.07 

2.0e-l 

6.00 

5 

1.3e-2 

1.07 

2.6e-F0 

1.50 

2.9e-2 

0.86 

1.5e-l 

0.39 

6 

7.0e-3 

0.93 

2.1e-F0 

0.27 

3.7e-2 

-0.35 

2.5e-l 

-0.71 

7 

3.3e-3 

0.99 

8.6e-F0 

-2.10 

2.1e-2 

0.86 

l.le-3 

0.23 


errors and convergence orders with p =1 and p~^ = 


Table 4. Example (1) 

10^, using method (2.7)-(6.2) with stabilization parameters 7 = 10. 


The results in Table show optimal asymptotical convergence for the errors ej^ and e^. 
In addition we observe a slightly sub-optimal convergence for the error e“, approximately 
1.8. The error does not seem to converge. As in the previous test the error 
converges asymptotically to 1 , however for the rest of the errors we do not observe 
convergence. 
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Figure 3. Approximate solution Example |(1)| (left) case F = dfi and for case 
F = (right). 
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Appendix A. Technical lemmas 


In this section we prove two technical lemmas involving geometrical estimates on elements 
intersected by the interface F. We remind the reader that we assume that F is a simple curve 
with an arc-length parameterization X : [0, |F|) —>■ F, and assumption (l)-(3) in Section 2.1 


We assume that r is the radius of our tubular neighborhood given in Lemma Equivalently, 
for any 0 < r < r we have 


Tub{T) = {x : dist(x,F) < r}. 

Moreover, for any x G Tub{r) there exists a unique xr G F such that dist(x,xr) = dist(x,F) 
and X — xr is perpendicular to F at xr • 


Lemma A.l. Consider the F^ a segment of the curve F and I the straight segment connecting 
the two end points ofFg. Assume that \^\ < r, where r is the radius of the r-tubular neighborhood. 
Then, it holds 

\Ts\ <2\£\. 

Proof. Let xi, X 2 be the endpoints of the line segment i. For i = 1, 2, let W be the line that is 
normal to F at Xj. Let S be the infinite region enclosed by Ni and N 2 , then let S = Tub{r /2)r\S 
be the tubular section of F^. We see that the area of S is given by |5| = r|Fs|. Let Mi i = 1,2 
be the two lines that are parallel to ^ and distance r from i. Consider the trapezoid, R, enclosed 
by Ni,Mi, i = 1,2. We note that the area of R is given by \R\ = 2r\i\. The result will follow if 
we show that S G R since this will imply that rlF^I < 2r\l\. 

To this end, we first note that the line segment I C Tub{r/2) since the distance of any point 
in £ to F is less than \ll2\ which we are assuming is less than r/2. Let x G S', then we know 
there exists a unique point xr G F^ where the distance from x to xr is less than r/2 and the 
line, which we denote by L^, that passes through x and xr is perpendicular to F at xp. We 
know that intersects I and we call this point y. We also know that distance between y and 
Xr is less than r/2 since y G £ C Tub{r/2). Hence, we have shown that dist(y, x) < r. This, of 
course implies that the distance of x to the infinite line iext (which is the line that contains the 
line segment i) is at most r. In other words, x is in between the two lines Mi and M 2 . Since x 
was in the tubular section S it was in between the two lines Ni and N 2 and hence x belongs to 
the trapezoid R. □ 



Figure 4. Illustration of a curve seg¬ 
ment Fg, straight line segment i con¬ 
necting the end points of F*, the tubu¬ 
lar section of radius S, and the trape¬ 
zoid R. 
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A.l. Proof of Lemma Let e = Int(9ri n 5 T 2 ) G with Ti,T2 G T^, and the previous 
definitions of = e H . We analyze the case in Q~. The same analysis is valid in 0+. We 
proceed analyzing two cases depending on the intersection of the edge e with T: 

(i) If e n r = 0, (e = e ). If either Ti or T 2 does not belong to the results is trivial. 
Consider the case where Ti and T2 belong to , the interface sections Ti^r and T 2 ^r are 
nonempty. We consider the midpoint me the segment e and the ball i?|e|/ 2 (Rie) of radius 
|e|/2 and centered at me- If the interface does not cross the ball B\e\/2{'<^e) then we have 

|Tf I > \B\e\/2{me) r^Tl\. 

Denote by a the minimum angle of the triangulation (given by shape regularity), and 
let d = min{a, 7r/4}. Now consider the isosceles triangle T) with base edge e and base 
angles a. Then, clearly Ti C i3|e|/2(jRe) C Ti and |Ti| = (|e“|/2)^ tanja}. Therefore, 
|7"rl > (|e“|/2)^tanjd}. Assume then that the interface crosses the ball B\^e\/ 2 {'f^e) 
in Ti- Observe that this implies that there exists a point in Ti^r who’s normal passes 
through me with distance less than |e|/2. Now, if r 2 ^r crosses the ball then we will have 
two points on r at a distance less than |e|/2 whose normal passes through me which 
contradicts the tubular neighborhood assumption. Therefore 

max \T~\ > f—'j tanldj. 

i={i,2}' * ' - V 2 / ^ ^ 

(ii) If e n r 7 ^ 0. Similarly as in the previous case, consider me the midpoint of e~ and the 
ball of radius |e“|/2 centered at me- We observe that this ball can not cross both Ti^r 
and T 2 ^r- Therefore 

/ I I \ 2 

max iTGl > (J-1) tanjdj. 

*={1,2}' * ' - V 2 / ^ ^ 


A.2. Proof of Lemma |ll[ We analyze the case in D . We first observe that, by triangle 
inequality we have 


\It\ 

i=l 


< 3 max |e 
*={1,2,3} 


Therefore, using Lemma A.l 


iTrl < 6 max |e,- 
*=1,2,3 * 


Same analysis is valid to prove the statement in D"*". 


Appendix B. Construction of basis functions of S^{T) 

In this section we define basis functions {rci, rc2,1^3} of the space S^{T), for T G . Consider 
{xi,X 2 ,X 3 } nodes of the element T. Let {Ai,A 2 ,A 3 } be the barycentric coordinates of a point 
X gT with respect to {xi,X 2 ,X 3 }- Consider the following representation 


Wi{x) 


Wi {x), 
wf{x), 


if X G T ; 
if X G r+. 


3 

i=i 
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We construct {wi,W2,W3}, a set of basis functions of the space S^{T), by satisfying the following 
conditions: for i = 1 , 2,3 


Wi{xj) = Sij 
[tCi(a:o)] = 0 

[DtoWi] = 0 

[DrigWi — 0 ] 


for j = 1 , 2 , 3 ; 
xq G Tr; 
to = t{xo); 
no = n{xo); 


where 6 ij 


1 , if i = j 
0, if i / j. 


These conditions are written in a system of size 6x6 using the barycentric coordinates repre¬ 
sentation of wf, i.e., we find the unknowns coefficients of Wi for i = 1 , 2 , 3 , solutions 

of; 


/ 

^l.+ 

0 

0 

^1.- 

0 

0 \ 


/ <1 \ 


f ki 

\ 


0 


0 

0 

^2.- 

0 


at,2 


Si ,2 



0 

0 

<^3,-1- 

0 

0 

^3.- 


a^3 


Si ,3 



Ai(a;o) 

A2(a:o) 

A 3 (a::o) 

—Ai(xo) 

—A2(a;o) 

~A 3 (a:o) 




0 



1 

b 

o 

—Rto A2 

— Dto^3 

Rto Ai 

DtgX2 

DtoXa 




0 


V 


A2 

P'^ DyIq A3 

0 

1 

P Rno A2 

p DjigX^ J 

V %3 / 


1 0 

/ 


where 



1, if Xj G 
0 , if Xj G T^. 
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